Site dilution in Kitaev's honeycomb model 
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We study the physical consequences ol site dilution in Kitaev's honeycomb model, in both its 
gapped and gapless phases. We show that a vacancy binds a flux of the emergent Z2 gauge field 
and induces a local moment. In the gapped phase this moment is free while in the gapless phase 
the susceptibility has the dependence x(h) ~ hi(l//i) on field strength h. Vacancy moments have 
interactions that depend on their separation, their relative sublattice, and the phase of the model. 
Strikingly, in the gapless phase, two nearby vacancies on the same sublattice have a parametrically 
larger x(h) ~ (/ifba(l//*)] ) 1 . In the gapped phase, even a finite density of randomly distributed 
vacancies remains tractable, via a mapping to a bipartite random hopping problem. This leads 
to a strong disorder form of the low-energy thermodynamics, with a Dyson- type singularity in the 
density of states for excitations. 

PACS numbers: 75.10.Jm, 75.40.Cx, 75.50,Mm 



I. INTRODUCTION 

The interplay of disorder and interactions is one of the 
most fascinating aspects of condensed matter physics. It 
is also one of the most challenging ones, and opportuni- 
ties for an exact analysis are rare. We show in this paper 
that Kitaev's honeycomb modeli offers a new and fruitful 
setting for such investigations. 

The introduction of randomness places the Kitaev 
model in the broader context of spin systems with 
quenched disorder, where the physics of spin glasses and 
infinite randomness phases are just two instance of con- 
ceptually new physics which have arisen. More specifi- 
cally, there is considerable interest in the question what 
happens when quenched disorder is introduced into a 
magnetic system exhibiting a novel correlated ground 
state, such as a classical or quantum spin liquid. Be- 
sides a search for new types of disorder-induced phases, 
it turns out that the properties of many quantum sys- 
tems can be probed sensitively through the controlled 
introduction of impurities. In particular, impurities may 
reveal elusive features of the clean system, as illustrated 
by the use of non-magnetic ions to uncover the order- 
parameter symmetry in a superconductor—. Experimen- 
tally, local probes such as nuclear magnetic resonance or 
atomic force microscopy can be used to distinguish im- 
purity from bulk susceptibilities?. 

An example of this behaviour in one dimension is site 
dilution of an antiferromagnetic spin-^ Heisenberg chain, 
which creates free chain end&2., leading to a Curie con- 
tribution to the susceptibility^. In higher dimensions, 
putative examples of spin liquids have in particular at- 
tracted attention, as for these no local diagnostic, such 
as an order parameter, is available for determining their 
nature. For instance, numerical work on perhaps the 
most enigmatic S — 1/2 Heisenberg magnet in d = 2, 
that on the kagome lattice, indicates that non-magnetic 
impurities generate a local dimcrisation pattern but do 
not induce a local moment. In an opposite extreme, the 



classical (gapless) spin liquid on the SCGO lattice ex- 
hibits a local moment evidencing classical fractionalisa- 
tion: in the low-temperature limit, its size is exactly half 
that of a free spin's^. This goes along with an extended 
spin texture visible in a modulated local susceptibility in 
NMR experiments^ a feature also predicted to occur for 
a candidate gapless quantum spin liquids on the kagome 
lattice^. 

On the level of a theoretical description, we are limited 
by the small number of instances of spin Hamiltonians for 
which a simple and controlled derivation of the existence 
and nature of a quantum spin liquid phase is available. 
Happily, the Kitaev honeycomb model provides a rare 
example of a solvable spin model with both gapped and 
gapless liquid phases^. As a particular attraction, solv- 
ability in this context implies not only detailed knowledge 
of the respective ground states, but also the availability 
of much information on the excited states. 

The Kitaev honeycomb model has attracted much in- 
terest in its own right as it exhibits some of the neces- 
sary elements to develop a quantum computer: it sup- 
ports fractionalized excitations, both Abelian and non- 
Abelian. There have been several proposals for an exper- 
imental realisation of the model, both with cold atoms in 
an optical lattice^ and in a solid state syste m 10 ' 11 . 

The solvability of the model has its origin in the ex- 
istence of an extensive set of non-dynamical fluxes of an 
emergent Z2 gauge field that permits a reduction of the 
Hamiltonian to a free fermion problem. We show in this 
paper that these steps remain tractable in the presence 
of site dilution and the Hamiltonian remains solvable, 
allowing us to calculate the magnetic properties of the 
Kitaev honeycomb model in the presence of vacancies. 

As one central result, we find that a vacancy binds 
a Z2 flux. Since these fluxes are the aforementioned 
anyonsi, this may have consequences for quantum com- 
putation: computations are performed by braiding the 
fluxes and if those braids encircle an impurity an addi- 
tional Aharonov-Bohm phase may be picked up. 
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Beyond this, gapped and gaplcss phases differ greatly 
in their response to the introduction of non-magnetic im- 
purities. In the gapped phase, we find that a single va- 
cancy generates a paramagnetic moment with a magni- 
tude that tends to that of a free spin as the gap becomes 
large. It is localised on one site adjacent to the vacancy. 
In the gapless phase, interactions with bulk excitations 
lead to an effective field-dependent moment, with a di- 
vergent susceptibility x(h) ~ ln( 1 //i.) , localised now on 
all three sites adjacent to the vacancy. 

The moments of different vacancies interact, in ways 
that depend on the phase of the system and the rela- 
tive sublattice of the vacancies. In the gapless phase, the 
most dramatic consequence arises for two nearby vacan- 
cies on the same sublattice. Here the impurity suscepti- 
bility is parametrically larger than for a single vacancy 
x(h) ~ (/i[ln(l//i)] 3 / 2 ) -1 . In the gapped phase, the in- 
teraction between vacancies on different sublattices de- 
creases exponentially with their separation. In a system 
with many vacancies we obtain a situation akin to the 
picture of the Bhatt-Lee singlet phase. This it turns out 
can be analysed as a random bipartite hopping prob- 
lem. We characterise the resulting broad distributions of 
energy levels underpinning the thermodynamics, with a 
density of states p(E) ~ F(E)/E, where T[E) vanishes 
slower than any power of energy E as E — > 0. 

In the remainder of this paper, we flesh out these asser- 
tions. We begin by introducing the model and its emer- 
gent degrees of freedom along with useful notation in 
Sec. HU Sec. IIIII presents the low-energy Hamiltonian in 
presence of vacancies and magnetic field. The gapped 
phase is discussed in Sec. lIVI A perturbative demonstra- 
tion of flux-binding is followed by successive treatment 
of the properties of the one-, two-, and many- vacancy 
problem. Sec. [V] is devoted to the gapless phase. Here 
we provide details for the analysis of the behaviour of an 
isolated vacancy and a vacancy pair. Sec. I VII summarises 
with a brief discussion and pointers to open questions. 
Two appendices treat some more technical material: how 
to project from the enlarged Majorana fcrmion Hilbert 
space down to the physical 5 = 1/2 Hilbert space; and 
details of the honeycomb lattice Green functions, which 
we use extensively. 

A short account of some of this work, particularly that 
on the gapless phase, has appeared in a Letter—, which 
also covered weak bond disorder. Some details of calcu- 
lations omitted from the present paper are described in 
Reffill Before proceeding, we alert the reader to a su- 
perficially different but conceptually related study of the 
coupling of a magnetic impurity to a spin in the Kitaev 
honeycomb model^. 



II. THE KITAEV HONEYCOMB MODEL 

The Kitaev honeycomb model consists of spins on the 
sites of a honeycomb lattice, with nearest neighbour in- 
teractions that depend on the orientation of the bond, 



unit cell 




(a) (b) 

FIG. 1: (a) Labelling (x,y,z) of the three bond types in a 
honeycomb lattice, and site numbering 1 to 6 for plaquette 
operator Wq. (b) The unit cell, lattice vectors and sublattice 
convention used throughout this work. 

labelled x, y or z as shown in Fig. [UJa). Representing 
a spin variable at site j by the Pauli matrices cr" and 
denoting the neighbouring sites by k, the Hamiltonian 
reads 

H = ~ £ J*5fd%- £ Jyttiil- £ J^5%. (1) 

x-links y-links z-links 

In the following we indicate the bond orientation between 
sites j and k using oijk = x,y or z. Without loss of 
generality, we take J a > 0. 

A. Mapping to free fermions 

This Hamiltonian is exactly solvable, for example with 
a local transformation cr" = ib^Cj which represents each 
spin using the four Majorana fermionsi and 
Cj. This transformation enlarges the Hilbert space and 
to emphasize this difference we mark operators in the 
Hilbert space of spins with a tilde and those in the space 
of Majoranas without. After this transformation the 
Hamiltonian reads 

Hu = g £ . ",!,<■:(■!, u jk = ib" ik III ■ . (2) 

The operators ujk commute with each other and with Hf t . 
One can therefore fix the values of (ujk) = Ujk = ±1, 
move to a subspace of the full Hamiltonian and obtain a 
bilinear form in the Cj 's. For each set {ujk}, the resulting 
Hamiltonian H u inherits a bipartite structure from the 
honeycomb lattice. This is displayed by introducing, for 
a lattice of N unit cells as depicted in Fig. [ljb), two 
A r -component vectors of Majorana fermion operators, ca 
and Cb, from sublattices A and B respectively, and an 
N x N matrix M, with entries J ajk Ujk- Then 
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The energy levels of H u can be expressed in terms of the 
eigenvalues ±5 m of the 2N x 2N matrix 
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The matrix H can be interpreted as a tight binding model 
and we will use this viewpoint extensively. The connec- 
tion between the eigenvalues of H and those of H u is as 
follows. The matrix M has a singular value decomposi- 
tion M = USV T , where U and V are N x N orthogonal 
matrices and S is an N x TV positive semidefinite diagonal 
matrix. Label a unit cell by r and let the two Majorana 
fermion operators within it be ca,t and ce. r . Then Ma- 
jorana fermions in the eigenbasis are 

Cm, A = ^ ' ^mr C A,r 
r 

and c„ ljB = ^V£ T c B ,r- (5) 

r 

Defining the complex fermions a m = -r ( .b) 
the Hamiltonian is brought into the diagonal form 

N N 

H u = i S m C m ,AC m: B = 2J ( 2a m a m ~ l) ■ (6) 
m— 1 m— 1 

The three Hamiltonians presented so far, H , and 
act in different Hilbcrt spaces and should not be 
confused. For instance, their respective Hilbcrt space 

dimensions are 2 2N , y/t^ N = 2 AN and 2 N . 
B. Emergent gauge field and non-dynamical fluxes 



operator must be applied to eigenstates of in order 
to obtain those of H. In fact, however, for the operators 
considered in this paper, matrix elements calculated us- 
ing eigenstates of H u are the same as those calculated 
using the projected physical states. A full discussion of 
the projection operation is given in Appendix [A1 

C. Excitations: fluxes and fermions 

The ground state energy in a particular flux sector is, 
from Eq. ©, E a = — J2 m ^mi an d the absolute ground 
state lies in the flux sector that minimises this energy. 
There are two distinct types of excitation. First, within 
the same flux sector, "matter" excitations involve the 
occupancy of the fermionic modes, (2d\ n a m — 1) = 1 
in Eq. <j6j> . Second, flux excitations consist of a set 
of flux values different from those in the ground state. 
For the disorder-free system, the ground state is flux 
freeJ^ As a function of the interaction parameters J Q , 
there are four phases in the absence of an applied field: 
three symmetry-equivalent phases in which matter exci- 
tations are gapped, when one exchange dominates (e.g. 
Jz > Jx + Jy), and one phase in which they are gapless, 
around the "isotropic" point J x =J y =J z =J . Flux exci- 
tations are generically gapped: their cost is given by the 
difference in fermionic "zero-point" energies, Eq, between 
excited and ground state flux sectors. 



Kitaev showed that there exist non- 
dynamical, commuting flux operators W = 
(Jj <J k <J k <7 l ctj <T m ""---CTn (T j denned along 
any closed loop on the lattice. They have eigenvalues ±1, 
and so are Z2 variables. Moreover, their effect on fermion 
hopping is that of a flux through the corresponding 
plaquctte. With periodic boundary conditions there are 
two independent non-contractible loops winding around 
the system. A complete and independent set of variables 
specifying the flux state is provided by the values of all 
but one of the fluxes through the hexagonal plaquettes, 
supplemented by those of the pair of fluxes through the 
non-contractible loops. 

Indeed, these fluxes are simply related to the variables 
Ujk encountered above. Numbering the sites around a 
plaquette from 1 to 6 [see Fig. [TJa)] the Z2 flux through 
a plaquette can be written as Wq = M21M23M43M45M65U61. 
The spectrum of H u depends only on these fluxes^ but 
because of invariancc under gauge transformations, many 
choices of the set {ujk} encode the same flux sector. The 
simplicity of the model is that the gauge field represented 
by Ujk has no dynamics. 

The appearance of a gauge degree of freedom is a con- 
sequence of the transformation to Majorana operators, 
which doubles the dimension of the Hilbert space for each 
spin. Using a variable 0j = ±1 on each lattice site, the 
gauge transformations are implemented in the standard 
way, as Cj — > OjCj and Ujk — > QjUjk9k- Due to the dou- 
bling in Hilbert space dimensions per spin, a projection 



III. INTRODUCING VACANCIES 

Vacancies have a dramatic effect on local properties 
of the Kitaev model. We will show that they bind a 
flux in the ground state. Furthermore, the removal of 
a site reduces the number of "constraints" on the spins 
adjacent to it and allows a magnetic moment to form. 
The formation of this moment, its susceptibility and the 
interaction between vacancy moments are the subject of 
this paper. 

Our analysis starts from the observation that the 
steps described above (representing spins using Majorana 
fermions and fixing the values of Ujk) remain valid in the 
presence of both vacancies and the leading terms in the 
Zeeman energy. The effect of a vacancy manifests itself 
on several levels in the description of the model. 

First, the three original plaquettes that meet at the 
vacancy site now form one big plaquette (different from 
all other plaquettes: see Fig. [2]) and hence the number 
of independent fluxes decreases by two. Second, and as 
a consequence, there is a qualitative change in the way 
that a Zeeman field couples to spin system. A key result 
we reach is a derivation of a modified version of the tight 
binding model, Eq. ([4]), valid in the presence of vacancies. 

To set the context it is useful to recall the effect of 
a Zeeman field in the model without vacancies.— On in- 
cluding the field, the Hamiltonian no longer commutes 
with the Z2 fluxes W, which therefore become dynamical. 
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FIG. 2: In the gapped phase (vertical bonds J z > J x + J y ) 
the ground state sector has a Z2 flux through the vacancy 
plaquette. 

This spoils solvability but a weak field may be treated 
perturbatively. Since there is a finite energy gap Aw 
between the ground state and other flux sectors, at field 
strength h -C one can project onto that flux sector. 
Practically all our results presented here thus pertain to 
a particular (typically the ground-state) flux sector. 

In more detail, for the undiluted lattice, matrix el- 
ements of the Zeeman energy between states from the 
same flux sector are all zero (see Ref. Q] and Appendix |A")) . 
The leading contribution to a projected Zeeman Hamil- 
tonian is therefore second order in h, which results in a 
non-vanishing but finite magnetic susceptibility at small 
hX 

The altered nature of the Zeeman coupling comes 
about because, with the merging of three plaquettes into 
one, four previously distinct flux sectors collapse into a 
single one. As a result, some terms in the Zeeman Hamil- 
tonian now commute with all fluxes W in the model with 
vacancies: the non-zero matrix elements of these terms 
connect eigenstates of H that belong to different flux sec- 
tors in the undiluted system but to the same sector in the 
diluted system. They thus generate contributions to the 
projected Zeeman Hamiltonian that are linear in h, and 
these are responsible for the dominant local magnetic re- 
sponse at weak field. Using the labelling indicated in 
Fig. [3|a), this part of the Zeeman energy is 

H z = -{Kdf + h y dl + h,5f) . (7) 

To show that the Hamiltonian H + Hz can be diag- 
onalised using the steps set out in Sec. [H] we find the 
counterpart to Eq. ([3]) arising from Hz- To start, recall 
that each bond of the lattice is associated with two Ma- 
jorana fcrmions. As removal of a site breaks three bonds, 
it leaves three unpaired Majorana fermions. We denote 
them by b", where j = 1, 2 or 3 labels the sites adjacent 
to the vacancy. With this notation, the contribution to 
Ha generated by the term had" in Hz is ihab'jCj. The 
Majorana fermions b" are represented in a tight bind- 
ing model by three orbitals that do not appear for the 
undiluted system. These orbitals are coupled to the ones 
representing the Cj 's by hopping of strength h a , as shown 
in Fig.^b). 




(a) (b) 

FIG. 3: (a) Sites at which field components h x , h y , and 
h z contribute linearly to the projected Hamiltonian with a 
vacancy, (b) Representation as an equivalent tight binding 
model, with three new sites coupled with respective hopping 
matrix elements h x ,h y ,h z . 



In studying this and related problems, it will be useful 
to keep in mind some basic properties of tight binding 
models with nearest hopping on bipartite lattices: with 
Na (Nb) orbitals on the A (B) sublattice, energy eigen- 
values generically consist of \Na — Nb\ zero modes and 
minjA^, Nb} pairs (the two members of a pair having 
energies of equal magnitude and opposite sign). 

The task then is to calculate the field dependence of the 
ground state energy of the Kitaev model with one or more 
vacancies, using the tight binding model of Fig.^b), and 
its generalisation in the case of more than one vacancy. 
We do this by using a T-matrix approach [see Eq. ([20])] to 
express properties of the system with vacancies in terms 
of the Green function of the hexagonal lattice tight bind- 
ing model, for which there are convenient analytic ex- 
pressions in the flux-free sector (reviewed in Appendix 
B). We will find that there are marked differences be- 
tween the gapped and gapless phases and we separate 
our analysis accordingly. 



IV. GAPPED PHASE 

We begin our discussion with the gapped phase, whose 
relatively simple structure allows us to make consider- 
able progress starting from the single vacancy problem. 
Indeed its very name holds the promise of a pertur- 
bative treatment, and for J z > J x + J y , we can use 
j x y = J x ,yl ' Jz as small parameters in an expansion. 

Our analysis proceeds in several steps. We first derive 
an effective Hamiltonian demonstrating that vacancies 
in the gapped phase bind a flux. Next, we provide a 
detailed analysis of the single vacancy problem, where a 
lone fermionic zero mode appears in the energy gap. This 
zero mode is localised in real space in a striking fashion: 
its probability density is zero outside a wedge that has 
the vacancy at its apex, and decays exponentially with 
distance from the apex. This mode we show carries an 
effective paramagnetic moment on the site linked to the 
vacancy by a strong bond; the size of the moment grows 
with decreasing j x .y 
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Zero modes belonging to vacancies on different sub- 
lattices "hybridise" when their wedges overlap, as we 
demonstrate by an analysis of the two- vacancy problem. 
For this, we can derive an effective Hamiltonian describ- 
ing the energetics in the gap in the presence of a field. 

Finally, the many vacancy problem leads us to consider 
an effective bipartite random hopping problem (BRH). 
This can be analysed in the spirit of a strong disorder 
renormalisation group treatment, as nearby pairs of va- 
cancies hybridise exponentially more strongly than dis- 
tant ones. This leads us to a strongly divergent low en- 
ergy density of states near E = 0. 

The remainder of this section gives a detailed account 
of this set of phenomena. 



A. Flux Binding 

The ground state flux sector of the undiluted Kitaev 
model is flux free, for example with Ujk = +1 on every 
link. However, as we now show, the removal of a site 
binds a flux to the vacancy plaqucttc. 

The energy differences between flux sectors may be 
found at small j x y without resorting to the Majorana 
decomposition given abovei Instead we extend to a sys- 
tem with vacancies the approach originally presented by 
Kitaev for the undiluted model. To this end, write the 
Hamiltonian as H = Hq + V, where 

H =-J z 22, °j°k 

z- links 

v = —j x E syss -J y E ■ 

x-links y-links 

The ground state of Hq is 2 JV -fold degenerate and has an 
energy Eq = —J Z N. The low energy states of H can be 
understood by projecting onto ground states of H and 
working pcrturbatively in V. In this subspace V takes 
the form of an effective flux Hamiltonian that acts within 
the ground states of H , lifting the degeneracy. 

We use standard perturbation theory to find the effec- 
tive flux Hamiltonian,— denoting the n th order pcrturba- 
tion by H^ s and with n a projection onto ground states 
of Hq. Further, let \a) and E a be the eigenstates and en- 
ergy levels of Hq. The action of V on such an eigenstate 
is to flip two spins and thus change the energy by 4J 2 . 
The first two terms in the perturbation theory are 

H$ = nvn = o 

£(2) _ V\a)(a\V u _NJ± (f + f) ^ 

a 

where the primed summation is over states outside the 
ground state manifold of Hq. At general order, H^ ^ 



is zero and the most important contribution to H^ s is 

' V\a)(a\V\b)(b\V---V\2n-l)(2n-l\V 
ab 2n-i (E -E a )(E -E b )---(E -E 2n -i) 

~2n (10) 

Here we omit for conciseness other terms in H^p' that 
have a subleading effect on degeneracies. Without dilu- 
tion, the lowest order term that reduces the ground state 
degeneracy is 

^=-^ 2 E^O+ const, (11) 

where Wq is the flux through a hexagonal plaquette of 
the lattice. The prefactor — is found from summing 
4! terms, corresponding to the permutations of elements 
from the four V operators. 

With dilution, at each vacancy there is one hexagonal 
plaquette that has a larger prefactor, indicated by • in 
Fig. [21 Here, in place of a spin pair with strong coupling 
J z , there is simply a lone spin. In consequence, when the 
two elements of V adjacent to the vacancy plaquette act, 
the energy change is only 2J Z . These elements are shown 
in Fig. [21 where we also define the vacancy plaquette G>. 

Perturbation theory must be extended to 8 th order to 
find a term that depends on the flux through the vacancy 
plaqucttes. This involves a summation of 8! permuta- 
tions from terms encircling £b and an additional 8! from 
those encircling both <£b and • plaquettes. The resulting 
Hamiltonian for an isolated vacancy reads 

+ j. tit E w 6 - E fa >< %) ) • 

(12) 

The fluxes Wq and Wm are determined by the fourth or- 
der terms in Eq. ([12]). Both the associated plaquettes 
are flux free in the ground state: (Wq) = (W+) = +1. In 
such states the terms involving Wq can be combined to 
give an energy contribution (9/2 11 ) J z jijy ^£v The 

positive coefficient of Y^, Wq indicates that the introduc- 
tion of vacancies changes the ground state flux sector: 
each vacancy binds a flux, so that in the ground state 
(W&) = — 1. This approach may be extended to include 
larger voids: the flux through those larger voids makes 
a contribution to the effective flux Hamiltonian that is 
higher order in j x and j y . Note that the resulting flux 
gaps are numerically rather small. For j x = j y = 1/3, 
one obtains a number of order 10 _6 J Z for a single va- 
cancy plaquette. The theory developed in the remainder 
of this section needs to be understood as applying for 
fields small on this scale. 
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B. A single vacancy 

Having seen that a vacancy binds a flux, we now dis- 
cuss its magnetic properties. For this purpose it is con- 
venient to switch techniques, from perturbation theory 
for the spin Hamiltonian to the Majorana fermion repre- 
sentation introduced in Section III Al We show that the 
vacancy generates a fermionic zero mode, which in turn 
leads to a twofold degeneracy of all eigenstates of H. The 
degeneracy is lifted by non-zero Zecman field h, reflect- 
ing the formation of a free moment distributed over sites 
close to the vacancy. The moment has an anisotropic 
effective g-factor that depends on the parameters j x ,y 

We first recall the gapped spectrum of matter exci- 
tations for the undiluted lattice, which can be found 
by Fourier transform of the tight-binding Hamiltonian, 
Eq. Q. We define basis vectors ni and 112 for the hon- 
eycomb lattice as in Fig. [TJa), and introduce a wavevec- 
tor q with components q\ = q.ni and q 2 = q.n 2 . 
The eigenstates of H with this wavevector have energies 
±5 q = ±| J z + J x e iqi + J y e iq2 \: these form two bands, 
arranged symmetrically around zero. The minimum of 
Sq over q in the phase under discussion (0 < j x ,jy < 1) 
is A = J z — J x — J y , yielding a band gap for H of 2 A. 

A single vacancy on one sublattice generates a zero 
mode, an eigenstate of the honeycomb lattice tight bind- 
ing model at zero energy with amplitude only on the op- 
posite sublattice. In the gapped phase the eigenfunction 
is exponentially localised and has a particularly simple 
form: sites with non-zero amplitude are located inside a 
60° wedge emerging from the vacancy and in the direc- 
tion of the strong bonds. The zero mode wavefunction 
has a straightforward expression in a system without flux. 
In this case, for a B vacancy at the origin, its amplitude 
on the A site in unit cell r is a representation of Pascal's 
triangle, with 

* B (r) = M(-l) n ^ j™ 1 j£ 2 ( ni n + r) ( 13 ) 

if the site lies within the wedge (ni,ri2 > 0). Thus, in- 
side the wedge the amplitude on A sites decays with dis- 
tance from the vacancy in a direction-dependent fashion. 
Outside this wedge, and on all B sites, the amplitude 
is zero. The wavefunction associated with a vacancy on 
the opposite sublattice is related to this one by inver- 
sion symmetry. Both cases are illustrated in Fig. 2J The 
normalisation constant 

n = \ (i ff 4f yiiii (14) 

will play in important role in what follows. In a system 
with a single flux bound to an isolated vacancy, the zero 
mode wavefunction has site amplitudes of the same mag- 
nitude as just described. The concomitant (Z2) phases 
are of course gauge-dependent. 

The presence of a zero mode is linked to the forma- 
tion of a free local moment around the vacancy, and this 
moment is polarised by the local field components that 
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A vacancy B vacancy 

FIG. 4: Vacancy zero modes from the tight binding model, 
Eq. Q. Area and colour of each circle respectively represent 
magnitude and sign of the wavefunction at that site. 

are included in projected Zeeman energy, Eq. (JTj)- To 
examine moment formation in detail, we compute the 
field-dependence of the ground state energy of the Kitaev 
model with a vacancy, using the tight binding model of 
Fig. [31b). As discussed, the basis orbitals for this tight 
binding Hamiltonian consist of all those appearing in the 
honeycomb lattice with a vacancy, together with three 
additional ones arising from unpaired Majorana fcrmions 

At h = 0, the matrix elements of the tight binding 
Hamiltonian involving these additional orbitals vanish. 
Its spectrum in this case therefore includes four zero 
modes located in the middle of the gap between positive 
and negative energy bands. The zero-mode subspace is 
spanned by the orbitals r x , r y , r z and the wavefunction 
#(r) (see Fig. 13(b) and Eq. (p]l ). 

At leading order h acts within this subspace, lifting the 
degeneracy of the zero modes. In fact, since the vacancy 
mode ^(r) has no amplitude on sites ri and r2, the or- 
bitals r x and r y are unaffected, ultimately yielding states 
with energies quadratic in h x and h v via coupling to the 
finite energy bands. 

By contrast, states arising from ^(r) and r z have en- 
ergies linear in h z for small h. Since these states lie 
within the energy gap, it is natural to project onto the 
two-dimensional subspace that they span. The projected 
tight binding Hamiltonian has the form 

\Nh z J' (15) 

where M appears as the amplitude ^(rs) on the site ad- 
jacent to the vacancy. Viewing the projection in terms 
of Majorana fermions (taking for definiteness the case of 
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a B vacancy at the origin) we have [following Eq. J5])] 

c M = ^ B (r)cA,r and c hB = K a . (16) 

r 

The Hamiltonian H u [Eq. ©] projected onto the low 
energy states can then be written as 

H u = 5i(2oJoi - 1) , (17) 

with Si = A/"|/i z |. The ground state magnetisation 
m z = —di lz E is thus J\fsgn(h z ) at leading order in h. 
Moreover, H u has a low-lying excited state, with the 
opposite magnetisation, which is higher in energy by 
2Af\h z \, formed by taking the occupation number a\a\ 
to be one rather than zero. 

We conclude that the vacancy has generated a para- 
magnetic moment of size J\f. A consequence of the form 
of the projected Zeeman energy Hi is that the magneti- 
sation associated with the moment is entirely localized 
on the site r3 adjacent to the vacancy. Indeed, for zero 
jx , jy ; the Kitaev model with a vacancy consists a number 
of spin pairs strongly coupled by J z exchange, together 
with one unpaired spin on this site. However, even in 
the limit of j x , j y small, the local moment is strikingly 
different from a free spin because the limits of small j 
and small h do not commute. In particular, irrespective 
of field orientation, only the z-componcnt of the moment 
develops a finite expectation value at small h, because 
all matrix elements of and er^ 3 within the zero-mode 
subspacc vanish. Equivalcntly, the local moment has an 
anisotropic g-tensor, with g zz as the only non-zero com- 
ponent, varying with j x and j y between g zz = on the 
phase boundary with the gapless phase, to g zz = 1 for 

jx,jy -> 0. 

This non-trivial form of the g-tensor reflects the fact 
that the local moment describes a collective coordinate. 
Some further insight comes from considering a system 
in which, additionally, the spin at has been removed, 
leaving vacancies on two sites adjacent in the z-direction. 
Without Zeeman terms, eigenstates of this spin model 
can be shown (we omit details) to have a two-fold degen- 
eracy that arises from the double vacancy, and we use |+) 
and |— ) to denote two chosen othonormal ground state 
wavefunctions. 

We employ the states |±) to construct the two lowest- 
energy eigenfunctions, |0) and |1), of single vacancy 
model including a weak projected Zeeman field. Intro- 
ducing the eigenstates | f) and | J.) of a^ 3 , these may be 
written as 

|0) = V(l + fe/2)| t) ® |+> + yffj. - fc/2)| I) ® |-) 
and 

|1) = v/(l - g zz /2)\ t) ® |-} + v/FW2)| I) ® |+) 

Then (0|er*j0) = g zz and (1|<t* 3 |1) = -g zz , but for any 
state in this subspace (ojL) = (c^ ) = 0. 



C. Two vacancies 

We next consider a pair of vacancies. Our main in- 
terest is in the interaction between the local moments 
formed near well-separated vacancies, and in particular 
we assume a sufficiently large separation that they are 
associated with different hexagons of the lattice. Then 
they generate two separate vacancy plaqucttcs of the type 
shown in Fig. [3] In the Majorana fcrmion representation 
of the Kitaev Hamiltonian, Eq. ([S]), two vacancies lead 
to six uncoupled Majorana fermions from the six bro- 
ken bonds. In the tight binding model, Eq. ([4]), these 
are represented as six uncoupled zero energy orbitals. 
At large vacancy separation, however, the leading Zee- 
man coupling involves only the z-components of field for 
jx ,jy < 1 j and so only two of these orbitals play an im- 
portant roloii. There are also two localised modes similar 
to that of Eq. (fT5|) . and so we must consider a total of 
four states in the low-energy subspace. 

The resulting behaviour depends crucially on whether 
the vacancies belong to the same or opposite sublattices. 
For vacancies on the same sublattice, the localised modes 
are both at zero energy. In this case, therefore, two 
vacancies produce two essentially independent paramag- 
netic moments, each similar to that for a single vacancy. 

By contrast, vacancies on opposite sublattices may give 
rise to moments that interact. More specifically, consider 
two such vacancies, placed so that each one lies inside the 
zero-mode wedge of the other. In this case the vacancy 
modes of the tight binding model hybridise, forming a 
pair of eigenstates at energies ±e, with e <C A if the 
vacancy separation is large. As in our discussion of a 
single vacancy, we include contributions to the Zeeman 
energy that are diagonal in the flux sector, with local 
fields hi and I12 acting near the two vacancies; when 
vacancy separation is large, we find again that only the 
z-components enter the energy at leading order, and we 
denote these by h\ and hi. 

We will show, extending Eq. (fT7|) . that H u for the two- 
vacancy problem at weak field, projected onto low energy 
states, has the form 

H u = Nhi{2a\ai-1)+Nh 2 (2ala 2 -1) 

+ ie(a\ + a 1 )(al + a 2 ). (18) 

As for a single vacancy, in the gapped phase with < 
jx,jy < 1) only the z-components of local moments 
develop non-zero values in weak fields. We can write 
spin operators at sites adjacent to each of the vacan- 
cies within a flux sector in terms of complex fermions. 
Projecting onto the four-dimensional space of low energy 
states, and denoting the projection operator by Q, we 
find Qaf n Q = N{2a} n a m — 1), with m = 1 or 2 labeling 
the vacancies. 

We next derive the effective Hamiltonian (fT8|) and cal- 
culate the energy e. We do this by using a T-matrix 
approach to relate the Green function for the tight bind- 
ing model that represents the system with vacancies to 
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the Green function for the undiluted hexagonal lattice. 
To establish some notation, consider matrices H, Hq and 
V, related by H — Hq + V, and define the Green func- 
tions G = (z — H)~ l and G = (z - Hq)- 1 . Then the 
T-matrix is 



and 



v(i-G vy 



G — Go + GqTGo 



(19) 



(20) 



In general, if the spectrum of Hq has an energy gap and 
H has levels within this gap, T has poles in the complex 
z-plane at the locations of these levels. We determine e 
by finding these poles. 

For clarity it is convenient to relate the undiluted lat- 
tice to the system with vacancies in two steps: in the 
first step we eliminate site orbitals at the locations of 
the vacancies; in the second we couple three additional 
orbitals around each vacancy, as illustrated in Fig. [3jb) . 
To implement the first step we choose V to be a potential 
1/e acting at each of the vacancy sites, and take the limit 
e -> 0. 

Consider a B vacancy in the unit cell at the origin and 
an A vacancy in the unit cell at r, with r = nini + n,2ii2. 
Let G^P (ri, ra) be the Green function for a lattice with- 
out vacancies but with fluxes through one of the plaquc- 
ttes adjacent to each of the sites where vacancies will be 
introduced. The poles of T are at the values of z for 
which 

G^(r,r)Gf B (0,0) - Gp B (r, 0)Gp A (0, r) = 0. (21) 

We are unfortunately able to compute Gp(ri, r 2 ) only 
numerically. We find however (in the sense made precise 
below) that key features of its behaviour are the same as 
for the lattice without fluxes. In the zero flux sector con- 
venient analytical expressions are available for the Green 
function, which we denote by Gq (0, r), and we base our 
initial discussion on these. We will see that e decreases 
exponentially with vacancy separation. To study well- 
separated vacancies we therefore require the Green func- 
tion at z small compared to the gap A. An expansion in 
powers of z/A gives (see Appendix IB")) 



G AA 



(0,0) 



j2 
Jx 



3y) 2 - ki'if 2 



f_1 \ni+n2 + l 

riBAfn \_ V tl Ani -n 2 (ni+n 2 \ 

U l U > r j- j Jx Jy { m ) 



o 



o 



V A ) 



n\,Ti2 > 0, otherwise 



(22) 

with the other Green function elements obtained by sym- 
metry. 

From this we find that two vacancies on opposite sub- 
lattices, located so that their zero-mode wedges overlap 



produce a pair of levels within the gap, with energies of 
magnitude 



m + n 2 
ni 



(23) 



when r is large compared to the decay length of 
the Green function. We also find, from a numerical 
study, that in the gapped phase Gp A (r, r) / 'G AA (r , r) and 
Gp A (0, r)/G AA (0, r) approach unity when |r| is large 
compared to the decay length of the Green function. Be- 
haviour of well-separated vacancy pairs is therefore the 
same in both the ground-state and the flux-free sectors. 

In a second step we include the basis orbitals in the 
tight binding model that arise from unpaired Majorana 
fcrmions 6q and b". As for a single vacancy, the leading 
contributions arise only from a = z when |r| is large. We 
are therefore concerned with four states within the gap, 
and provided Zeeman fields are sufficiently small we can 
project onto this subspace. Doing so, we obtain a tight 
binding model of the form 



M T ^) With M= ( A ? 1 Nh 2 
The Majorana fermion Hamiltonian 



H u = i ( b A c A )M 



cb 



(24) 



(25) 



can readily be rewritten in the form of Eq. (fTS)) . It has 
eigenvalues ±S+ and ±SL given by 



S 4 



N /£ 2 +AA 2 (/i 1 + / l2 ) 2 - y /e 2 +Af 2 (h 1 -h 2 )' 



and so Eq. (|T8|) has the form 



H u 



S + (2a\ a+ - 1) + 5_(2a f _a_ - 1) 



(26) 



For example, with h\ = h 2 = h, this gives a moment for 
the system MS M whole, of 



m(e, h) = 



4j\f 2 h 



y/e 2 + Ah 2 N 2 



(27) 



For fields /i<e there is a large, field independent impu- 
rity magnetic susceptibility \ = 4A/" 2 /e, and for e <C h 
the induced moment saturates at 2Af, twice the isolated 
vacancy moment. In Fig. [5] wc illustrate behaviour for 
hi y= h 2 . This figure also demonstrates that projection 
onto states within the gap provides a very accurate treat- 
ment of the full Hamiltonian with Zeeman energy of the 
form given in Eq. ([7]). As for a single vacancy, from this 
form of the Zeeman energy it is apparent that the mo- 
ment is entirely localised on the sites equivalent to r 3 in 
Fig-OHb), adjacent to each vacancy. 
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bipartite hopping Hamiltonian 



-0.4 



h 1 x°100 



FIG. 5: Local moments formed near each vacancy, as a func- 
tion of field strength, in a system with two vacancies on op- 
posite sublattices: m z vs. hi for (from the top) /12 = e/2 and 
h,2 = e/6. Interaction strengths j x = j y = 1/3; vacancy sepa- 
ration 10(111+112); mode splitting e = 7.9 x 10 -3 J z ; saturation 
moment M = 0.863. Points are obtained from the projected 
Hamiltonian, Eq. (|18[) ; curves are for the full Hamiltonian in 
the flux-free sector using lattice of size 20 x 20. 



D. Finite density of vacancies 

We next investigate the properties of the Kitaev model 
in the gapped phase with a finite density of vacancies. We 
will see that couplings between the low-energy degrees 
of freedom generate an impurity band with a density of 
states that has a Dyson-type divergence at zero energy. 
This leads to a divergent macroscopic susceptibility. 

The approach to the two vacancy problem outlined 
above can be extended to include a finite density of va- 
cancies. Consider a system with Na vacancies on the 
A-sublattice and Nb vacancies on the B-sublattice, and 
let Gp A , G BB and G AB be matrices of Green function 
elements for the undiluted lattice (all functions of z) eval- 
uated between the vacancy sites. The energy levels of the 
honeycomb lattice tight binding model with these vacan- 
cies are given by the values of z for which 



G AA (G BA ) T 



C<BA 



= 0. 



(28) 



We assume that vacancy separations are much larger 
than the decay length of the Green function. Then 
Na + Ng vacancies give rise to an impurity band of 
Na + Nb levels, with energy width <C A, at the centre of 
the band gap of the undiluted lattice. In these circum- 
stances we are concerned with \z\ -C A and the Green 
function elements can be approximated by the leading 
order in an expansion in z/ A. At this order Gp B is eval- 
uated at z = and 



/~<AA 
Kj p 



z[d z G F (0,0) 



1 



-(z/JV 2 J*)-l, 



and similarly for G BB . Within these approximations, the 
impurity band levels are therefore the eigenvalues of the 





(~<BA 



(G BA 




(29) 



Our focus is on the case of compensated vacancies (Na = 
Nb) for which the impurity band is generically free of 
zero modes; by contrast uncompensated vacancies result 
in at least \Na — Nb\ zero modes. 

The Hamiltonian of Eq. ([2"9")) is one example from the 
well-studied class of bipartite random hopping (BRH) 
models. Such models are characterised by a density of 
states p(E) that is strongly divergent as energy E ap- 
proaches zero, having the form p(E) = J 7 (E)/\E\, where 
T(E) is a function that goes to zero more slowly than 
any power of E but ensures that the density of states is 
intcgrablei 18 ' 19 While renormalisation group treatments 
of these models typically generate broadly distributed 
coupling strengths, even our bare Hamiltonian already 
has couplings that vary over many orders of magnitude 
since they depend exponentially on vacancy separation. 

We find from a numerical study of Eq. (|29|) a strong di- 
vergence in the density of states, over about ten decades, 
as illustrated in Fig. [Sj While the exact behaviour of 
F(E) is difficult to ascertain, our results are compatible 
with the form 



T(E) cx 



and a fit yields x = 1.7. 



1 



l0g[l/£]* 



(30) 




-2.5 
-ln[ln[1/E]] 

FIG. 6: Density of states for 160 randomly placed compen- 
sated vacancies on a lattice of 3200 sites, with fluxes attached 
to vacancies and with J x = J y = 1, J z = 4. The density of 
states is an average over 1000 disorder realisations. Solid line 
is a best fit to the data with gradient of 1.7. Inset: The same 
data, plotted on a scale that illustrates the very large energy 
range considered. 

We next examine response to a Zeeman field. Con- 
sider first the impurity band Hamiltonian in the absence 
of a field. The off-diagonal blocks of -ffeRH have a sin- 
gular value decomposition G BA = vsu T , where (taking 
Na = Nb) v and u are Na x Na orthogonal matrices and 
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s is an Na x Na diagonal matrix with positive entries Sj . 
The impurity band Hamiltonian is therefore reduced to 
a direct sum of 2 x 2 block-diagonal forms by the trans- 
formation W t HbhhW, with 



W 



v 
u 



(31) 



Coupling to the projected Zeeman held [Eq. ([7|)] is in- 
variant under this transformation, provided the held 
strengths h a acting near each vacancy are the same. For 
this reason, each pair of levels from the impurity band 
with energies ±Sj makes a contribution to the Zeeman 
response like that from a single pair of vacancies with 
e — Sj. Following Eq. (|27|) . the total magnetisation and 
susceptibility are then 



X 



m(E,h)p(E)dE , 

4M 2 E 
(E 2 +4A/" 2 /i 2 ) 3 / 2 



T(E)dE . 



(32) 



The susceptibility integral has its largest contribution 
from a region around \f2Mh. Since F{E) is slowly vary- 
ing, the susceptibility has the form 



X 



(33) 



This singular behaviour arises because well-separated va- 
cancy pairs give rise to local moments that are fully po- 
larised even in a weak held. 



V. GAPLESS PHASE 

The response of the Kitaev model with vacancies to 
a Zeeman field is very different in the gapless phase 
compared to that in the gapped phase. The difference 
arises because the impurity modes that form a finite- 
dimensional low-energy subspace in the gapped phase be- 
come continuum resonances in the gapless phase. As a 
physical consequence, magnetic response in the gapless 
phase has striking singularities at weak field, which we 
discuss in this section. 

Behaviour is qualitatively the same throughout the 
phase, and we focus on the isotropic point J x =J y =J z =J. 
We will be concerned with the magnetisation of a sys- 
tem with a single vacancy in a Zeeman held that is weak 
compared to J, and with properties of a pair of vacan- 
cies that have a separation large compared to the lattice 
spacing. A feature of the gapless phases is that the zero 
mode induced by a single vacancy without Zeeman cou- 
pling is merely power-law localised in the gapless phase, 
and indeed is not normalisable in an infinite system.— 
Moreover, its probability density is not confined to a 
wedge (as in the gapped phase), but instead approxi- 
mately isotropic. One direct consequence is that there is 
a large response to all components of the Zeeman field, 
rather than just h z as is the case for < j x , j y < 1. 



We first comment on the ground state flux sector in a 
system with a vacancy. The approach used above (see 
Section IIV Ap to discuss flux binding to a vacancy deep 
in the gapped phase is not useful in the gapless phase 
because the expansion parameters j x ,jy are not small. 
Instead, as we have reported elsewhere^ a direct nu- 
merical calculation can be used to show that a vacancy 
binds a flux with energy — 0.027J. 



A. Magnetisation of a single vacancy 

In this subsection we calculate the magnetic response 
of the Kitaev model with a single vacancy in the gapless 
phase. We present results for both the ground state flux 
sector, with a flux bound to the vacancy, and for the 
flux-free sector. The second case has the advantage that 
calculations are simpler. In addition, it turns out to be 
relevant to behaviour of a system with two vacancies that 
both have bound fluxes: since we are dealing with Z2 
fluxes, low-energy properties in a system with two nearby 
fluxes arc like those in the zero flux sector. 

A summary of the main results is as follows. In the 
ground state flux sector the vacancy susceptibility at 
weak held h and low temperature T diverges as 



X <* 



\n{l/h) for 
ln(l/T) for 



T< h 
h<£T. 



(34) 



In the zero flux sector, the vacancy susceptibility has the 
still stronger singularities 



X oc 



l/[/i(ln(l//i)) 3 / 2 ] for 
l/[Tln(l/T)] for 



T4Z h/y/\n (l/h) 



(35) 

To derive these results, we require the ground state 
energy or (at finite T) free energy of the fcrmionic degrees 
of freedom in the relevant flux sector. These follow from 
the eigenvalues of the tight-binding Hamiltonian, and the 
necessary information is contained in the trace of the 
Green function and its dependence on the complex energy 
z. We express this Green function in terms of the one for 
a system without vacancies using the T-matrix approach 
described above [see Eq. (19)}. 

Let G(h,r,r') be the matrix element between sites r 
and r' of the Green function for the tight binding model 
of Fig. [3£b). We dehne the difference between the Green 
function trace at finite field h and at zero field, as 



p(z, h) = Tr[G(h, r, r') - G(0, r, r')] 



(36) 



The discontinuity in the imaginary part of p(z, h) across 
the real z axis gives in the usual way the difference be- 
tween the hnite and zero held density of states for the 
tight binding model, and so the ground state energy dif- 
ference is expressed by the integral 



2,771 



zp{z, h)dz , 



(37) 
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taken in an anti-clockwise direction on a contour enclos- 
ing the negative real axis. Similarly, the free energy dif- 
ference, computed within the flux sector, is 



T(h) = — — I Tp{z, h) In [2 cosh(z/T)] dz . 
2m J 



(38) 



evaluated on the same contour. 

We next outline the evaluation of £(h), taking h = 
(0, 0, h) for simplicity of presentation. In this case we can 
omit the orbitals arising from the Majorana fermions b x 
and b y , retaining only the one from b z (see Fig. [3]). We 
start with a lattice of 2N sites before dilution and give the 
& 2 -orbital the site label r z . To obtain the Green function 
G(0,r,r') for a system with a vacancy, we introduce a 
potential 1/e at the vacancy site r„ and take the limit 
e^O. In this limit G(0, r, r') consists of one (2N — 1) X 
(2N — 1) block, with elements 



G (r,r')- 



Go(r,r t ,)Go(r ^ ,,r , ) 
G (r v ,r v ) 



(39) 



for r, r' ^ r v ,r z and one lxl block with element 1/z 
for the site r 2 . We include non-zero h by applying the 
T-matrix approach a second time, with G(0, r, r') as the 
initial Green function, taking 



V 



h 
h 



(40) 



in the basis of sites r z , r 3 (see Fig.[3|). With the shorthand 
g(z) = G(0,r3,r3), the T-matrix is 



T 



h ( g(z)hz z 
z - g{z)h 2 \ z h 



(41) 



Finally, using J2 r G (0, r, r 3 )G(0, r 3 , r) = -d z g(z), we ob- 
tain 

p(z, h) = h 2 \z- x g{z) - d z g(z)]/[z - h 2 g(z)} . (42) 

Setting g(z) — zd z g(z) = a(z) + \b{z) and z — h 2 g(z) = 
u(z) +iv(z), with a(z), b(z), u(z) and v(z) real, Eq. (|37|) 
becomes 

8 {h )= h -f d,^tt Wx) . (43) 



u 2 (x) + v 2 (x) 

To make use of these results, we require matrix ele- 
ments of the Green function for the undiluted lattice. At 
small h, the function u{x) has a zero at x = xq, with 
— 1 <C xq < 0, and the dominant contribution to the in- 
tegral in Eq. (|43| comes from the vicinity of this point. 
We are therefore concerned with the Green function at 
small z: it has the behaviour (see Appendix IB"]) 

G (r,r)~ Azln[-( M z) 2 ] G (r 3 , r„) ~ -v , (44) 

where, for J a = 1, A = 1/\^3tt and /i = v = 1/3. The z- 
dependences of Eq. (|4^|) are a direct consequence of the 
masslcss Dirac spectrum for the nearest neighbour tight 



binding model on the honeycomb lattice, and hold with 
appropriate values for A, \x and v throughout the gapless 
phase. 

We expand the integrand of Eq. (|4"3")l about x = Xo in 
the small-/i limit, retaining only the leading terms. Using 
the asymptotic forms for Green function elements given 
in Eq. (gtj), we find 



u(x) « x + 



2,, 2 



K 2 v 



2\x In : 



and therefore 



hv 



X 



yj2\\xi{l/h) 



(45) 



(46) 



Then d x u(x)\ x=Xo ~ 2, and writing x = xq + s we have 
?lh\~ K2 r a a ( x oH x o) h 2 a(x ) 

^^L ds wR~' (47) 

Moreover, h 2 a(xo) ps 2xo, which yields the energy 



£{h) 



^2Aln(l//i) 



the magnetisation 



m(h) = —dh£{h) 



y/2\\n{l/h) 



(48) 



(49) 



and the susceptibility, Eq. (J35J). Results for non-zero tem- 
perature are obtained in a similar way. 

Behaviour for a general field orientation can be ob- 
tained from a similar, although more involved, calcula- 
tion. We find that, even in the general case, only the 
field magnitude enters the leading contribution to p(z, h) 
at small h, which is therefore orientation-independent. 
As for the gapped phase, the magnetisation is entirely 
localized on sites adjacent to the vacancy but now each 
of these sites, labelled ri,r 2 ,r 3 in Fig. [3Jb), carries 
a separate component (m x ,m y and m z , respectively), 
proportional to the corresponding component of h: as 
the field orientation changes, the induced magnetisation 
moves around the vacancy in real space! 

We next turn to behaviour in the ground state flux sec- 
tor, with a flux attached to the vacancy plaquette. All 
relevant information is contained in the function g{z), 
and in particular, its form for small z. This in turn 
depends on site-diagonal and nearest-neighbour Green 
function elements. Moreover, since the nearest-neighbour 
Green function elements are finite at z = in both flux 
sectors, the crucial quantity is the site-diagonal element, 
or equivalently the local density of states (LDOS). In the 
zero flux sector this varies linearly at small energy, re- 
flecting the Dirac cones in the honeycomb tight binding 
model dispersion. We use a numerical study to find the 
LDOS at small energy in the presence of a flux. The 
geometry employed and the results obtained are illus- 
trated in Fig. [7] We evaluate the LDOS at sites close to 
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FIG. 7: Results from a numerical study of the effect of a 
flux on the local density of states (LDOS). Bottom: geometry 
of the system, showing two fluxes (indicated by 7r's), joined 
by a string. The sign of the hopping matrix element on the 
links crossed by the string is reversed in the presence of the 
fluxes. Top Right: detailed view of the geometry, showing one 
flux (n) through a plaquette, and the labeling (a, b, c and d) 
of nearby sites. Top Left: LDOS as a function of energy E. 
Points: with a flux, at the four labeled sites as indicated; line: 
without flux (same behaviour at all sites). 



a plaquette threaded by an isolated flux, for a cylindrical 
system of infinite length and finite circumference (taken 
to be 1502 lattice units for the data shown). We find 
that the flux has a dramatic effect, resulting in a finite 
LDOS in its vicinity at small energy. In consequence, 
whereas g(z) is divergent for z — > in the zero flux sec- 
tor, it has a finite limit in the ground state sector. The 
resulting weak-field form for the susceptibility is shown 
in Eq. (|34| . 



B. Zeeman response with two vacancies 

The response to a Zeeman field of a system in the gap- 
less phase with a pair of vacancies involves the physics of 
the Kitaev model in particularly rich ways. In summary, 
we find three types of behaviour in the ground state flux 
sector, depending on the field strength h, the vacancy 
separation, and whether vacancies lie on the same or op- 
posite sublattices. The magnitude of the vacancy sepa- 
ration vector d sets a field scale h c , which decreases with 
increasing separation, (i) For fields large on this scale the 
response of a pair of vacancies is simply the sum of the 
responses that would arise for each in isolation. By con- 
trast, weak field behaviour depends on the relative sub- 
lattices of the vacancies, (ii) Two vacancies on opposite 
sublattices generate a contribution to the susceptibility 



that is field-independent at h <C h c and paramctrically 
larger at large |d| than the susceptibility per spin of the 
host system. (Hi) Most strikingly of all, two vacancies 
on the same sublattice have a susceptibility that is para- 
metrically more strongly divergent for h <C h c than for 
an isolated vacancy, being of the same form as for a sin- 
gle vacancy in the zero flux sector, as given in Eq. (|35p . 
In essence, this is because two vacancies bind a pair of 
Z 2 fluxes in the ground state sector, and because in their 
influence on low-energy properties these fluxes effectively 
fuse and cancel. 

We derive these results by expressing the energy 8(h) 
of the system with vacancies and a field in terms of Green 
function elements for the undiluted lattice, using exten- 
sions of the T-matrix methods outlined above. In the 
zero flux sector, analytical expressions for these elements 
are available, which we supplement with computational 
results in the ground state flux sector. Since calculations 
are quite involved, we omit many details and consider 
only the case of projected Zeeman fields [Eq. ©] ori- 
ented along the z-axis and of equal strength h for both 
vacancies. 

As in our discussion of a system with a single vacancy 
in the gapless phase, it is convenient to separate calcu- 
lations into two steps, introducing vacancies at the first 
step, and including a Zeeman field at the second. Some 
notation is summarised in Fig. [5] The Green function 
G(0, r, r') describing a system with N unit cells after the 
first step is a block diagonal matrix consisting of two 
lxl blocks with entries 1/z, from the 6| orbitals, and 
one (2N - 2) x (2N - 2) block with entries 



G(0,r,r') = G (r,r') 

+ (G (r,r vl )G Q (r,r v2 ))T v 



G (r vl ,r') 
G (r„ 2 ,r') 



for r,r' ^ r„i ,r v2 , r zl or r z2 , the form of the T-matrix 
being 



T„ = - 



G (r vl ,r vl ) G (r v i,r v2 ) 
G (r v2 ,r vl ) Go(r v2 ,r v2 ) 



(50) 



From this, using the definition given in Eq. (|36|) . we find 

p(z,h) = N(z,h)/D(z,h) (51) 

with (after a lengthy calculation) 

N(z,h) = 2/i 2 {z" 1 /i 2 [G(0,r 1 ,r 2 ) 2 

-G(0,n,ri) 2 ] + G(0,n,ri)} 

- 2/ i 2 [z-/ l 2 G(0,r 1 ,r 1 )]9 z G(0,r 1 ,r 1 ) 

- 2/i 4 G(0,r 1 ,r 2 )S z G(0,r 1 ,r 2 ) 

and 

D(z, h) = [z- h 2 G(0, n,n)] a - h 4 G(0, n, r 2 ) 2 

where we have used the symmetries of G(0, r, r') to sim- 
plify expressions. 
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FIG. 8: Site labelling used in our calculation of the Zeeman re- 
sponse of a system with two vacancies on opposite sublattices. 
Similar notation is used for vacancies on the same sublatticc, 
where ri and r2 are the sites adjacent to the vacancies in the 
z-direction, selected by our choice of field orientation. 



Behaviour in each of the cases summarised 
above can be extracted by considering simplifications of 
Eq. (f5Tj) in the relevant limits. First we note that the con- 
tributions to N(z,h) and D(z,h) involving G(0,ri,r 2 ) 
may be omitted when discussing weak field behaviour, 
because they appear only in terms that are 0(h A ). With 
this simplification, p(z, h) reduces to twice the expression 
that applies in a system with a single vacancy [as given 
in Eq. I[i2"|)] but with G(0, ri,ri) evaluated for the two 
vacancy system. We therefore need to examine the be- 
haviour of G(0, ri,ri) = g(z) in this case. The integrals 
in Eqns. (|37|) and (f38|) arc dominated by contributions 
from \z\ ~ 0(h) (omitting for simplicity logarithmic cor- 
rections), and so we are concerned with g(z) at the scale 
\z\ ~ h. We first discuss this and its consequences in the 
zero flux sector. 

(i) Independent vacancies. Coupling between vacan- 
cies in the expressions for G(0, r, r') involves Go(r„i, y V 2 ) ■ 
This is exponentially small in |d| unless \z\ < |d| _1 , 
which leads us to identify h c = J/\d\: for h S> h c we 
can neglect Gq(t v i,r v2 ) and the response is a sum of 
independent contributions arising from each of the two 
vacancies, (ii) Coupled vacancies on opposite sublattices. 
For h <C h c and compensated vacancies, we find that 
g(z) ~ |d| 2 z if \z\ <C |d| _1 and is small otherwise. In 
consequence, we obtain a field-independent susceptibil- 
ity that is of order |d| and hence much larger than the 
susceptibility per spin of the undiluted system. (Hi) Cou- 
pled vacancies on the same sublattice. In this case, using 
the forms for the Green function of the undiluted lattice, 
discussed in Appendix B, we find the same singularity 
in g(z) at small z as for a single vacancy, and therefore 
obtain the same singularity in the susceptibility. 

As a final step, wc consider g(z) in the ground state 
flux sector, with fluxes through each of the two vacancy 
plaquettes. Because the presence of fluxes breaks transla- 
tional invariance, we no longer have analytic expressions 
for the Green function of the undiluted lattice. Instead, 
wc calculate this, and hence g(z), using a numerical im- 
plementation of the T-matrix approach to obtain results 
for a pair of vacancies and fluxes, with separation in the 
range 100 - 200, embedded in an infinite lattice. We 
write <?Q^ = G(0, ri,ri), with a, ft labeling the vacancy 
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FIG. 9: Comparison of the dependence on complex energy z 
of the Green function element g^ B in a system with two va- 
cancies on opposite sublattices. Top: in the flux free sector. 
Bottom: in the ground state flux sector. Vacancies have sep- 
aration d in the x direction, with values of d as indicated. 
The asymptotic form in the flux free sector at small z is 
g£ B oc d 2 z\n(-z 2 ). 



sublattices. In Figs. [9] and [10] we compare the small- z 
behaviour of g^f in the ground state and zero flux sec- 
tors, for vacancies on opposite and the same sublattices, 
respectively. These results demonstrate that the form 
of the Green function at small z is indeed the same in 
both flux sectors, as anticipated from the idea that for 
low energy properties nearby fluxes effectively fuse and 
therefore cancel. The main distinction between the two 
flux sectors is a large scale factor relating behaviour in 
the two cases. 



VI. CONCLUSIONS AND OUTLOOK 

In summary, this work builds on our observation that 
the Kitaev honeycomb model is solvable in the presence 
of vacancies and the leading Zeeman coupling. Indeed, 
Kitaev's original solution strategy - to identify a non- 
dynamical flux field and analyse an effective hopping 
problem for fermionic variables in each sector - remains 
applicable. 

This is all the more remarkable as several observables 
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FIG. 10: As Fig. [9] but for vacancies on the same sublattice. 
The asymptotic form for small z in the flux free sector is 
g£ A oc [z\n(-z 2 )]-\ 



change in a fundamental way in the presence of vacan- 
cies. Firstly, the ground state ceases to be flux free: in- 
stead, vacancies each bind a unit of the emergent Zi flux. 
Secondly, the finite linear local susceptibility of the pure 
system is replaced by that evidencing the formation of a 
local moment around the vacancy, with an entirely free 
moment whose size varies continuously with the coupling 
constants in the gapped phase. Such moment formation 
can happen because a vacancy locally reduces the num- 
ber of constraints on the spins adjacent to it. Thirdly, 
the vacancy moments interact, which is somewhat in con- 
trast to the ultra-short range spin correlations of the pure 
system. In the gapped phase, we discover the physics of 
strong disorder fixed points via a mapping to a bipartite 
random hopping problem. Most remarkably, in the gap- 
less phase, the magnetic response of two nearby vacancies 
on the same sublattice is parametrically enhanced with 
respect to the single-vacancy case. The resulting forms 
of the susceptibility are not only a test of our assertion 
that vacancies bind a flux but also of the basic descrip- 
tion of the Kitacv model in terms of the variables laid 
out in section |TTJ 

In the absence of an experimental realisation of Ki- 
taev's honeycomb model, it might seem premature to 
ask how one would go about experimentally probing the 
phenomena we have described. Given the intense inter- 



est in realising this system via a cold-atom simulator—, 
this question is perhaps not quite so far-fetched, par- 
ticularly in view of recent developments towards single- 
site microscopy in optical lattices. These may open a 
unique window on introducing and locally probing the 
quantum state around vacancy degrees of freedom inter- 
acting through a strongly correlated bulk. The detection 
of the bound emergent flux presents a particularly 
exciting challenge. 

In addition, with a conventional condensed matter set- 
ting in mind, our results illustrate for this model how dis- 
order can serve as a probe of quantum correlated matter. 
Local probes of magnetism, such as NMR or muon spin 
rotation, have a parametrically enhanced local suscepti- 
bility around the vacancy compared to the bulk. More- 
over, since the local susceptibility of a vacancy depends 
on its flux state, these probes even give sensitivity to the 
flux degrees of freedom. 

This behaviour is fundamentally distinct from what we 
found in the case of slowly varying, weak bond disorder—. 
There, the main consequence of the presence of disorder 
was not in the magnetic response but rather in the heat 

2 

capacity, C oc T~>-+ & , which exhibits a downward drift 
in its exponent with disorder strength parametrised by 
A cx (SJ) 2 . Yet other types of disorder are compara- 
tively featureless - in the case of ± J bond disorder, the 
zero-field spectrum remains entirely unchanged as that 
disorder can be 'gauged away' by placing a flux on every 
plaquctte with a negative product of J's around it. 

Plenty of open issues remain. For instance, a descrip- 
tion of the many- vacancy behaviour based on an effective 
Hamiltonian would clearly be desirable, particularly in 
the gapless phase, as would be a simple formulation of 
the fermionic low-energy excitations of the system with 
a single vacancy and its attendant flux. 

Finally, there is as yet no commonly accepted classi- 
fication (or, indeed, definition) of quantum spin liquids. 
Whereas for the case of gapped spin liquids, a classifica- 
tion via emergent gauge fields at least seems reasonable, 
the situation is quite unresolved in the case of gapless 
ones. Indeed, one of our result presents a step backwards: 
since bond disorder leads to a continuous drift of expo- 
nents, it is not even possible to diagnose a 'Dirac spin 
liquid' via heat capacity exponents without knowledge 
about presence and nature of random strains. Nonethe- 
less, the richness of phenomena induced by static vacan- 
cies does suggest that in working towards a systematic 
understanding of spin liquids, the response to disorder 
of a magnetic state of interest may present particularly 
valuable clues. 
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Appendix A: Projection and the physical subspace 

In this appendix we discuss the projection operator 
that is necessary to obtain physical states of the sys- 
tem from eigcnstates in the enlarged Hilbcrt space in- 
troduced with the transformation from spins to Majo- 
rana fermions. In the extended Hilbert space of H u , 
Eq. (j5J), each state has at least a degeneracy of 2 2N , as 
can be seen in the following manner. Define an oper- 
ator D 



j 



+1 for each lattice site in the 
Hilbcrt space of spins. In the Hilbcrt space of Majoranas 
this operator is -D, = b^b^Mcj and commutes with the 

3 3 3 •* 

Hamiltonian H u There are 2N operators Dj, each 
with eigenvalues ±1, which leads to the aforementioned 
degeneracy. 

Knowing that D = +1, it is clear that physical states 
of the Majorana system must have positive eigenvalue 
for Dj on every site. Kitaev constructed a projection 
operator 



l + Dj 



(Al) 



that ensures this. 

In the remainder of this section we review an a naly sis of 
the projection operator, expanding on Refs. 2lH22| . and 
also show that for the observables we consider, matrix 
elements evaluated using eigenstates of the Hamiltonian 
H u are the same as those obtained using the projected 
physical states. To further our efforts in understanding 
the effects of the projection operator, it is instructive to 
form complex fermions from the Majoranas b x ,b y ,b z ,c: 
with sites j, k nearest neighbours on the A, B sublattices 
respectively, define new variable s 21 ' 22 



OLjk 

3 



ib. 



(x?) f = \ (&; 



■1 " 

ib. 



fr = 5 (CA,T + iCB,r) fl = \ [c A ,r - iCB,r) ■ 



The Xr are located on the bonds of the lattice and / r 
in the unit cells, as shown in Fig. [TTJ In this notation, 
let site j be in a unit cell at r. Then Ujf. = ib a ^ h b°^ k = 

2(x") t Xr - !• Eigenstates of Ha, Eq. ©, are direct 
products of a wavefunction \x) for the gauge degrees of 
freedom, and a wavefunction |/) for the matter fields: 
1$) = \x) ® |/}- The choice of \x) encodes the flux sector 
and gauge, while |/) is an eigenstate of the Hamiltonian 
H u , Eq. ((HI), which in the notation of Eq. (|A2|> takes the 
form 



M A = | (M T - M) M s = i (M T + M) 



H u = (P f) 



M A 
-M s 



M s 
-M A 



f 
f 



(A3) 



In terms of the complex fermions x an< i /; f ne opera- 
tors D j now take different forms on the A and B sublat- 
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FIG. 11: The re-fermionization of Majorana fermions to com- 
plex fermions. Variables Xr are located on the bonds of the 
lattice and f T in the unit cells. 



tices: 

D A . r = (X? f + X X r)(X V r f + Xl)(xl 1 + XlWl + fr) , 

D B , r = (Xr-l, - X?_„ 1 )(X^- t „ 2 - XiU 2 )(Xr Z f ~ X r )x 

iti-fr)- 

(A4) 

The projection operator can also be re-written as 



V 



1 

$2N 



/i,r<i/,r' 



fi.r 



(A5) 

The operator D A ^ r changes the bond fermion num- 
ber on the three bonds attached to site t A (b), leaving the 
flux sector unchanged. It also acts on the matter state 
|/). Notice, however, that acting with Dj on the sites 
at both ends of a bond leaves the fermion number on 
that bond unchanged, and acting with Dj on all sites in 
the lattice leaves all bond fermion numbers unchanged. 
Defining D = FJ . Dj and V as the sum of all operators 
in V that change the bond fermion number in an inequiv- 
alent way, normalised by 1/2 2JV_1 since there are 2 2Ar_1 
terms in V ', the projection operator can be rewritten 2 ^ 
as V = V'(l + D)/2. Here D gives the parity of the total 
fermion number: with iV v and Nf the number of bond 



(A2) 

and matter fermions respectively 



D = (-lf*(-l) 



(A6) 



In this form, it is clear that the projection operator anni- 
hilates states of odd total fermion number. The complex 
fermion Hamiltonian H u in Eq. (|A3|) conserves fermion 
number modulo 2 and therefore commutes with D, which 
also commutes H u . Contrast this with Dj that commutes 
with H u but not with H u . One can then block diago- 
nalise H u into blocks that act on Hilbert spaces of even 
and odd fermion number. 

To summarise, the projection operator either annihi- 
lates a state that has odd total fermion number or trans- 
forms a state with even total fermion number to an equal 
weight superposition of all terms in V that change bond 
fermion number in an inequivalent way- 
Having established the effects of the projection op- 
erator we will now show that for a large class of op- 
erators, matrix elements evaluated using an eigenstate 
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states. 



Appendix B: Green's functions of the hexagonal 
lattice 



1. Gapless phase 



FIG. 12: A graphical illustration of the final state that results 
from the action of various spin operators on an initial state 
\ X ) with all (u jk ) = (2 ( X ?y Xr - 1) = +1. We denote u jk = 
+1 by black lines and Ujk = —1 by thick red lines. Spin 
operators act on both the gauge and matter fields but we 
illustrate here only the gauge fields, (a) A single spin operator 
er™ changes the bond fermion number on a single bond and 
changes the flux sector, adding it flux either side of the bond, 
(b) and (c) Nearest neighbour spin operators a" erf change the 
bond fermion number and flux sector unless they are of the 
form cr . Jh o k Jk . (d) A three spin operator that changes bond 
fermion number but not the flux sector and can by undone 
by a term in V' ■ Care must be taken with operators of this 
type as the terms in V' also act on the matter sector |/). (e) 
At the edge of the lattice, spin operators do not necessarily 
change the flux sector. 



1$) = |x) ® I/) °f Ha are the same as those obtained 
using the projected physical states. 

It is immediately apparent that if an operator O 
changes the fermion numbers on the bonds in a manner 
that cannot be undone by a term in V' then its expec- 
tation value is zero. This is the case for a single spin 
operator or for a two spin operator that is not nearest 
neighbour and in the direction of the bond: if O = c^er^ 

then (O) = in the ground state unless j, k are nearest 
neighbours and /? = 7 = ajk- The Kitaev honeycomb 
model therefore has only nearest neighbour spin corre- 
lations non-zero.— The effects of various spin operators 
are illustrated graphically in Fig. [T21 

For simplicity it is desirable to calculate matrix ele- 
ments using an cigenstatc |5>) = |x) <8> |/}- It is only per- 
missible to do so if these matrix elements are the same 
as those obtained using the projected physical states. 
Throughout this paper we consider operators O that 
leave unchanged the bond fermion number. This includes 
the Hamiltonian and types (c) and (e) in Fig. [T2J For 
this class 

($\v6r\<i>) (<$>\dv\<s>) ($|6(i + £>)|$) 



($|7>P|$) 



(&\v\&) 



($1(1 + £>)|*) 



($1$) 

(A7) 



At the steps of this derivation we have used sequentially 
the following facts: all spin operators commute with V 
and V 2 = V\ only the identity part of V leaves the bond 
fermion number unchanged; and |$) is an eigenstate of D. 
For operators of this type we are thus free to evaluate ma- 
trix elements using an unprojected eigenstate and obtain 
the same result as when we use the projected physical 



We seek the Green function for the hexagonal lattice at 
small energies. With uniform exchange coupling J a = 1, 
transforming Eq. Q to momentum space, the Hamilto- 
nian and corresponding Green function are 



Go(q) 



/q 
/q 

(2 = 



Z /q 
/q/q" V /q Z 



(Bl) 



where / q = (l+e lqi +e tq2 ). The real space Green function 
is then 

Gf (r) = J G^(q)e-^0. (B2) 



If z is outside of the band then Gq (r) can be evaluated 
directly for certain elements and analytically continued to 
the whole of the complex plane.— Here r denotes the unit 
cell and the sublattice indices a, /3 are given explicitly. 
Using complex energy z = E — ie, with E real and e a 
positive infinitesimal, we define 



.4 



y/(z-mz + 3) 



D 



16z 



(z-l)3(z + 3) 



(B3) 



It can be shown that22 



zAK 



-zA(K-2iK') 



E > 1 



E < 1 , 



(B4) 



where K = K(B 2 ) is the complete elliptic integral of 
the first kind, with complex parameter k 2 and K' = 
K(l — k 2 ). For small z we find the asymptotic form 
given in the text. Using the symmetry of the lattice and 
the defining Green functions equations, it can be shown 
that zG£ A (0) - 3G^ B (0) = 1. The Green function at 
an arbitrary site may be found from recursion relations 
involving sites closer to r = 0^ 

We now discuss the Green function at sites far from the 
origin. For z« 1 and r» 1 the dominant contribution 
to the integral (|B2[) is from the region of small |/ q |. We 
use the fact that the spectrum |/ q | is asymptotically lin- 



ear close to the Dirac nodes ±Q, where Q 



2 71 



2 71 
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Let q = ±Q + p, then 



j2 _ | p 2 



2 4tt 2 



?2 _ 3_2 



3z /•Im[(p :c +ip y ) e - i Q- r ] _ 4p . r d 2 



,2 _ | p 2 



2zi z2rz . 
= ]Bm (Q.r- 



4tt 2 
(B5) 



where 8 is the angle between r and the x axis and K\ is 
first order modified Besscl function of the second kind. 
Similarly, it can be shown that 



G^(r) ]sin(Q.r 



(B6) 



_Jf [_]cos(Q.r), (B7) 



where i*To is the zeroth order modified Besscl function of 
the second kind. 



z small compared to the gap J z ~J x —Jy can be found by a 
perturbative expansion in j x and j y . Let r = 711111+712112 
and / q = (J 2 + J x e i?1 + J^ 2 ). Then 



- l/qP 



4tt 2 



= -/jT 1 (l+i.e-^+ive-^)- 1 e-*»*0 

= (_l)|ni| + |n 2 |+lj-l j-lml j>2| f I" 4 ' + l 712 '^ 

2 x y V M / 



n ii n 2 < 0, otherwise. 



V ^1 



ni,n2 > 0, otherwise. 



(B8) 



(B9) 



2. Gapped phase 

We consider the parameter regime J z > J x + J y > 0. The Green function between sites on the same sublattice 
The Green function for the gapped phase at real energies can be obtained in an analogous manner. We find 



G^(r)^(-l)l ni H" 2 l +1 — jKI x 
for sgn(nl) = sgn(n2) 

for sgn(nl) 7^ sgn(n2) 



E 

fei,fe 2 =0 



E 

fci,fc 2 =0 



ni| + \n 2 \ +h + k 2 \ (k x + k 2 \ , 2kl . 2k2 
\n 1 \+k 1 )\ h I Jx Jy 



nil +h + k 2 \ /|n 2 | +ki +k 2 \ 2ki , 2fe2 

|m| + *i /V M + Ai h 
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